Simulation of Muon Tomography Projections to Image the Pyramids of Giza

Purpose: A geometric simulation of a possible two-plane detector was developed to test the abilities of the detector to generate high-resolution images of the Great Pyramid using muon tomography. Methods and Materials: Trajectory range, angular resolution, and acceptance of the detector were calculated with a simulation. Trajectories and the corresponding sinogram space covered were simulated first with one detector in one location, and then two moving detectors on adjacent sides of the pyramid. The resolution at the center slice of the pyramid was calculated using the angular resolution of the detector. Results: The simulation returned trajectory range encompassing the pyramid and peak angular resolution of .0004sr. Sinogram space covered by one position was inadequate, however two moving detectors on adjacent sides of the pyramid cover a significant portion. Resolution at the center of the pyramid is roughly 3m. Conclusions: The simulation provides a way to calculate the detector positions needed to cover an adequate amount of sinogram space for high-resolution cosmic-ray tomographic reconstruction of the Great Pyramids.


Introduction
Knowledge of the Great Pyramid of Giza is currently limited to knowledge from optical instruments.The question of how these structures differ, and reasons why, call for detailed, high-resolution imaging of those internal core structures by other means.In 2017 a new void was detected by the Scan Pyramids team through cosmic-ray muon tomography, albeit with relatively poor resolution [1].A project has since been proposed by an interdisciplinary group of researchers from the Oriental Institute, the University of Virginia, Fermilab, the Geophysics Department of Cairo University, and the University of Chicago Medical Physics collaborating to produce high-resolution tomographic images.These would reconstruct the interior structure of the Great Pyramid of Giza with one-sided cosmic-ray muon tomography [2].The Medical Physics aspect involves exploring the ability of the proposal to meet the required tomographic sampling necessary for a high-resolution 3D reconstruction of the pyramid.This rotation project focused on building a simulation to study the geometric dimensions of the proposed detector and the Great Pyramid, angular resolution and acceptance of the detector, resolution of detection, and the range of sinogram space able to be covered.

Methods and Materials
The dimensions of the proposed detectors and the dimensions of the external shape of the Great Pyramid were taken from a previous GEANT4 monte carlo simulation of cosmic-ray muon tomography [3].A numerical simulation was written to show that the widest trajectory angles subtended the entire Great Pyramid if the far detector side was located 25m from the edge of the pyramid.Given the pixel dimensions of the detector, the angular resolution and acceptance of the detectors were calculated using the definitions provided by Lesparre et.al. [4] From this, a simulation was written to translate all potential trajectories of the detector to sinogram space covered to test the ability of high-resolution tomographic reconstruction.This simulation can be used to determine placement and movement of the detectors at the base of the pyramid.Lastly, the simulation was extended to know the resolution at the furthest depth of the pyramid given the angular resolution of the pixels.As this is a purely geometric simulation, the cosmic ray spectrum, probability of interaction, energy of incident muon, scattering, density, and attenuation were not included.These sections are explored in detail in the subsections below. 2

Geometric Dimensions
The proposed detector is built within a shipping container, of length 9.6m, height 2.4m, and depth 2m [3].A front panel and a back panel of scintillator modules allow hits to be connected to define a muon's trajectory through the detector.On each panel 2 × 2cm 2 pixels are arranged such that there are 120 pixels vertically, and 480 pixels horizontally.This creates a grid of 240 × 480 = 115200 ij-pixels (i representing the horizontal location, j representing the vertical location), each of which have a trajectory through the matching 115200 ij-pixels D = 2m away on the opposite panel.This generates 1.3 × 10 10 paths between the two panels.Trajectories of muons are calculated by reconstructing the vector connecting the pixel in the front array and back array knowing the distance between the two pixels m = ∆i and n = ∆j with ⃗ v = [∆i, ∆j] as shown in Figure1.The trajectories with the greatest angle provide the limit of the angular range.The pyramid with dimensions of depth and width 230.33m and height 138.7m must fit within this angular range to be imaged in its entirety.

Angular Resolution and Acceptance
Given the dimensions and number of the detector pixels, the angular resolution and acceptance can be calculated.Using the definition from Lesparre et.al. [4] the angular resolution for each discrete direction of sight of the detector is equal to the solid angle covering all trajectories able to hit a given pair of pixels of that direction of sight.This angular resolution was simulated with and without tilt correction.This angular resolution and the area of the detector with this angular resolution were used to calculate the acceptance.
The direction of sight r mn is the direction of the vector passing through the center of two pixels.For pairs of pixels with the same horizontal and vertical displacement (m = ∆i, n = ∆j), the direction of sight is equal.The angular resolution for this particular direction of sight is given by the solid angle between these two pixels δΩ = A/r 2 = 4cm 2 /r 2 , as shown in Figure 2. Correcting for tilt, for each direction of sight r mn the angular resolution is given by the solid angle δΩ = A cos θ cos Ψ r 2 . This takes into account how the area covered by the solid angle formed by the connection of the far corners of offset pixels (m = ∆i ̸ = 0, n = ∆j ̸ = 0 is not equal to the area of the pixels, as shown in Figure 3. With the solid angle as a function of the discrete directions of sight available from the detector, the number of pixel pairs that have that discrete direction of sight can be converted to the detection area available per solid angle.This is done with Equation 23 from Lesparre et.al. [4], τ (r mn ) = S(r mn ) × δΩ(r mn ) representing Acceptance = Detection Area × Angular Resolution for each discrete direction of sight r mn with m = ∆i, n = ∆j.Doing this for every r mn returns the acceptance τ (r mn ) of the detector.

Sinogram Space
Looking at one level of pixels, meaning all pixels of a row with a constant j (∆j = 0), the trajectories measurable through a 2d slice of an xy layer of the pyramid at the height corresponding to j can be calculated.The pyramid coordinate system is constructed with the origin at the center, x parallel to the i-axis, y parallel to the depth of the detector, and z normal to the surface of the earth.From this single horizontal slice the angle of a trajectory is given by ψ = tan − 1(dy/dx) where dy is equal to the 2m for all pixels due to the distance between the two detector plates, and dx = ∆i.Converting from detector ij space to pyramid xyz space we can form projection space ξ and η with ξ as the axis rotated by the azimuthal angle ϕ = ψ − π 2 [5] around the z axis.In projection space the angle of rotation ϕ as well as the intersection of the trajectory with the ξ axis allows calculation of that trajectory's location in a sinogram.A plot of projection space of one trajectory superimposed on the pyramid and pyramid coordinate system is shown in Figure 4.
With the left bottom corner of a panel being the origin of ij space, and the angle of rotation of η being ψ = ϕ + π/2, ξ = x cos ψ + y sin ψ and η = −x sin ψ + y cos ψ with ψ being the angle of the muon trajectory, i.e. the projection.The origin is shifted from the detector coordinate basis to the pyramid's coordinate basis with x = (i − (960/2))/100 and y = −230.33/2+ 25 both now in units of meters.With these conversions, for every possible muon trajectory representing a different angle and position of a projection, the sinogram space covered by the detector can be calculated.A plot of a subset of projection angles is shown in Figure 5 as example of the range of projection angles for one row of the detector at one position.
Finally, as the amount of sinogram space covered by the detector in one position can be calculated, the detector can be moved and rotated in pyramid space and the amount of sinogram space covered by different positions, and multiple detectors, can be shown.As the pyramid is stationary the time of projection images captured can be ignored.

Resolution
Lastly, the solid angle associated with each pixel pair was projected forward to its resolution at the center of the pyramid.The distance travelled and the resulting resolution is dependent on which vertical slice of the pyramid is imaged, with the center being the worst resolution assuming 360 degree rotation of the detector about the pyramid.For simplicity the lowest row of pixels was chosen, j = 0, as it has the most relevant angle range, and that row was paired with every pixel on the front detector.This was done for one position of the detector with significant overlap of solid angles for different pairs of pixels with the same direction of sight.

Geometric Dimensions
The geometric simulation of the maximum range of a stationary detector placed at x = 0 and 25m from the base of the pyramid contained trajectories that encompassed the furthest edges of the pyramid.This is shown in Figure 6.The normalized direction of sight for every eighth pixel on the far panel with every eighth pixel on the front panel is shown in Figure 7 with the detector facing up.

Angular Resolution and Acceptance
The angular resolution, i.e. solid angle, associated with a certain direction of sight, is plotted as a function of θ and Ψ without and with tilt correction as shown in Figures 8 and 9.These values are as expected, as Lesparre et.al. finds a peak angular resolution of .015srwith 5x5cm 2 pixels and a D = 80cm, because (5×5) (80/2) 2 ) = .015sr[4].Similarly, with the proposed detector with 2 × 2cm 2 pixels and a D = 200cm, the peak angular resolution is (2×2)  (200/2) 2 ) = .0004sr,which agrees with the simulation.Acceptance, i.e. detector area times angular resolution for each direction of sight, is shown in Figure 10

Sinogram Space
The points of sinogram space covered by the detector at one position are shown in Figure 11.This is reasonable as the width of the detector is 9.6m, while the width of the pyramid is 230m.Therefore at ϕ = 0, completely forward direction of sight, the total range of ξ 0 covered is equal to the width of the detector.When a detector is moved in 20m increments along the x-axis of the pyramid at a constant 25m from the base, and another detector is moved along the y-axis of the pyramid, a majority of sinogram space is covered as shown in Figure 12.

Resolution
Lastly, the projections and the resulting angular resolution are shown as a grid at the center slice of the pyramid in Figure 13.There is a resolution of roughly 3 × 3m 2 which leads to an oversampling at every 2cm pixel.

Discussion
This simulation showed the medical physics side of muon tomography and the possibility of high-resolution reconstruction using muon trajectories as tomographic line integrals.It also studied the angular resolution and acceptance of a current version of the detector, which can be changed based on resolution demands.The geometric simulation can also be used to test possible movements of detectors around the pyramid and the resulting sinogram space covered.The simulation shows that with detectors moving along two sides of the pyramid, a majority of sinogram space can be covered allowing possible reconstruction.The sampling resolution of sinogram space wanted for high-resolution reconstruction determines the step sizes between detector placement.The resolution at the center of the pyramid is roughly 3 × 3m 2 , leading to oversampling.
The main limitation to the simulation is that it is purely geometric, does not take into account possible scattering, energy attenuation, the muon spectrum, timing, or cosmic ray statistics.Future work involving statistical reconstruction involves taking into account coulomb scattering and scatter based on material properties as discussed by Schultz et.al [6] and Hanson et.al. [7].Specifically the energy loss incurred by muons of initial energy E 0 > E min travelling through the pyramid can be used to determine the opacity along the light of sight.That is, Hanson et. al.[7].
This means that given the path length L, muons of a certain incident energy E 0 lose an amount of energy per density, and that knowing that path length one can know the opacity the muon travelled through to yield a final energy at the detector of E f .The muon spectrum measured by the detector of open air could return a good estimate of E 0 , and the change in the resulting spectrum per muon trajectory through the pyramid can be used to calculate the opacity of that trajectory.
This simulation provides a way to calculate the positions needed to cover an adequate amount of sinogram space for high-resolution tomographic reconstruction, as well as a way to calculate the path length L through the pyramid as a function of trajectory.In future work the energy attenuation as a function of opacity, comparison to open air, and scattering should be included.

Conclusion
This simulation of geometric tracks of muon trajectories through a pyramid to a detector used in one-sided muon tomography showed the ability to fill sinogram space with a series of positioned detectors along the basis of the pyramid.It also showed the angular resolution and acceptance of the detector, and the resolution at the center of the pyramid.This suggests that using muon attenuation to calculate opacity of every trajectory one can reconstruct a highresolution image of the pyramids using one sided tomography.
. The values calculated in Lespare et.al. are not matched by our detector.The peak acceptance calculated by hand for a 16 × 16 pixel grid with pixel size 25cm 2 with peak angular resolution of .015sr is equal to 16 × 5 2 × .015= 100.However the figure presented by Lesparre et.al. has a peak of 200cm 2 sr.The plot shown in Figure 10 follows the calculation done above, rather than the plot shown in the Lesparre paper.

Figure 1 :
Figure 1: A diagram of the shipping container detector proposed using the ij coordinate system.A trajectory and its vector is shown.

Figure 2 :
Figure 2: A diagram of the un-tilted solid angle and calculation used in the code.

Figure 3 :
Figure 3: A diagram of the tilted solid angle and the calculation used in the code.

Figure 4 :
Figure 4: A diagram of projection space ηξ from one incident muon with trajectory ⃗ v, angle ψ, traveling through the bottom slice of the pyramid shown in black, with the pyramid coordinate system xy for reference.

Figure 5 :
Figure 5: A diagram of the projection angles captures by every sixteenth pixel in a single row paired with every sixteenth pixel in the front panel of the same row.If every pixel is used the 'fan' is filled in.

Figure 6 :
Figure 6: A plot of the range of trajectories with respect to the pyramid.If placed at the middle and 25m from the base of the pyramid, the entire pyramid falls within the trajectory range.

Figure 7 :
Figure 7: A plot of the directions of sight as unit vectors using every eighth pixel.

Figure 8 :
Figure 8: A plot of the angular resolution of the detector in sr per angle of direction of sight.The effect of tilt on the solid angle is not taken into account.

Figure 9 :
Figure 9: A plot of the angular resolution of the detector in sr per angle of direction of sight.The effect of tilt on the solid angle is taken into account.

Figure 10 :
Figure 10: A plot of the acceptance of the detector in cm 2 sr per angle of direction of sight .The effect of tilt on the solid angle is taken into account.

Figure 11 :
Figure 11: A plot of the points in sinogram space covered by the detector located at the center of the x-axis and 25m from the base of the pyramid.

Figure 12 :
Figure 12: A plot of the points in sinogram space covered by a detector moved along the entire range of the pyramid along both the x and y axis eleven times ( 20m steps).In both cases the detector is pointed perpendicular to the axis it is travelling along, meaning it is facing the pyramid.

Figure 13 :
Figure 13: A plot of the resolution at the center slice of the pyramid for every eighth pixel for reference.